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Abstract 

In this paper, we reexamine the possibihty of shear-driven turbulence during sed- 
imentation of dust grains in the protoplanetary disk. The shear-driven turbulence is 
expected to occur before the onset of the gravitational instability for MMSN model. 
While according to previous studies without taking account of growth of dust grains, 
with the larger abundance of dust grains, the gravitational instability is indicated 
to occur before shear-driven turbulence. In this paper, the case with dust growth is 
considered, and it is found that the Kelvin-Helmholtz instability tends to occur be- 
fore the gravitational instability even in the case with large abundance of dust grains. 
This is different from previous results without the dust growth. 

Key words: instabilities — methods: numerical — planetary systems: proto- 
planetary disk — solar system: formation 



1. Introduction 

For the formation of planetesimals, growth and motion of dust grains in the protoplane- 
tary disks are essential. In a typical scenario, dust grains settle toward the midplane, and a thin 
dust layer is supposed to form (Nakagawa et al. 1981; Nakagawa et al. 1986). The gravitational 
instability (GI) of a disk (Safronov 1969; Goldreich & Ward 1973; Sekiya 1983) is expected to 
occur when density in the disk exceeds the critical density that is given by 
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I3 ■ 



(1) 

where r is the heliocentric distance (Sekiya 1983, 1998). In order to resolve the problem of 
radial drift of meter-sized dust (Adachi et al. 1976), GI in the thin dust layer is essential. 
However, as dust grains settle toward the midplane, vertical dust density gradient increases. 
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In such a case, the vertical shear of the rotational velocity in the dust layer becomes strong. 
The strong shear is expected to induce the Kelvin-Helmholtz instability (KHI) (Chandrasekhar 
1961). As a result, KHI is expected to induce shear-driven turbulence. It is considered that 
turbulence due to KHI prevents dust grains from settling further toward the midplane (Sekiya 
1998; Sekiya & Ishitsu 2001; Michikoshi & Inutsuka 2006). Sekiya (1998) shows that GI does 
not occur in MMSN (minimum mass solar nebula) model (Hayashi 1981; Hayashi et al. 1985). 
On the other hand, Sekiya (1998) and Sekiya & Ishitsu (2001) show that turbulence induced 
by KHI becomes weaker if dust abundance in the protoplanetary disk is much larger than 
MMSN model. If turbulence is weak, GI will occur and planetesimals may form. To discuss 
the possibility of the shear-driven turbulence, Richardson number is known as an indicator of 
KHI (Chandrasekhar 1961). If the Richardson number is smaller than a critical value, KHI is 
expected to occur. Chandrasekhar (1961) showed that the critical value is 0.25, but Gomez & 
Ostriker (2005) and Johansen et al. (2006) showed that the inclusion of the Coriolis force yields 
a much higher critical Richardson numbers. Other than KHI, the dynamics of the midplane has 
been shown to be dominated by the streaming instability (Youdin & Goodman 2005; Johansen 
& Youdin 2007). 

Investigations of the turbulence are essential for understanding the processes of the 
planetesimal formation, and there are many previous studies (Weidenschilling 1980; Cuzzi et al. 
1993; Dobrovolskis 1999; Johansen & Youdin 2007; Bai & Stone 2010). These previous studies 
confirm the result that shows that GI is expected to occur before KHI if dust abundance is 
large. However, above previous studies did not take account of dust growth. Not only gas 
turbulence but also dust growth are essential to understand the planetesimal formation. Dust 
grains grow due to dust-dust collisions while they settle toward the midplane (Nakagawa et al. 
1981; Nakagawa et al. 1986). Nakagawa et al. (1986) shows that the law of gas drag force on 
dust grains is changed in the process of sedimentation of dust grains because of dust growth. 
They also show that the change of the law of gas drag influences sedimentation of dust grains. 
Dust growth is strongly related to gas drag. However, it is difficult to simulate both dust 
growth and gas turbulence numerically because of the computer performance. 

In this paper, we numerically calculate the sedimentation and the growth of dust grains, 
and the possibility of KHI is discussed using the distribution of dust density that is consistent 
with their sedimentation in the disk. In §2, the models of gas and dust grains of the protoplan- 
etary disk in this paper are summarized. In §3, we show numerical results for sedimentation of 
dust grains without growth. In §4, results are shown for the case with dust growth. In §5, we 
discuss the effects neglected in this paper. In §6, we summarize our results. 
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2. Models of the protoplanetary disk 



2.1. Disk model 

For simplicity, we restrict ourselves to the case with r = 1 AU. The protoplanetary disk 
is composed of gas and dust. In this paper, the gas surface density Eg and the dust surface 
density are assumed to be 

Sg = 1.7xl0Vg(Y[Xu]) (2) 

Sd = 7.1/d6ce(Yj^) '[gcm-2], (3) 

where /g, /d and .^ice are parameters of abundance for gas, dust and condensed water ice, 
respectively (Hayashi 1981; Hayashi et al. 1985). At r = 1 AU, .^icc = 1. The case with /g = 1 
and /d = 1 corresponds to MMSN model. In this paper, we consider the case with /g = 1 and 
various /d. In MMSN model, temperature T is given by 

2.2. Gas properties 

In this paper, we assume that the gas component is in hydrostatic equilibrium in the 
vertical direction without self-gravity. In this case, the gas density is given by 



V^^g 
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(5) 



where z is the height from the midplane of the disk and is the scale height of the disk given 
by 

^..^.4.7X10-(^)'|AU1. (6) 

The symbol fix is Keplerian angular velocity. The symbol Cg is the sound velocity as 



where k-Q is the Boltzmann constant and m^{= 3.9 x 10^^^ g) is the mass of gas molecular (the 
mean molecular weight is 2.34). The mean free path of gas molecules, Zg, is given by 1^ = 1.44 cm 
at r = 1 AU (Nakagawa et al. 1986). From equation (6), the protoplanetary disk is geometrically 
thin. 

2. 3. Dust properties 

In this paper, dust grains are assumed to be compact and spherical with radius s. In 
the case with s < 3/g/2 = 2.2 cm, the gas drag force is given by Epstein's law. In this case, at 



r = 1 AU and z = 0, the stopping time tstop is given by 

where ps is the internal density of the dust grain. In the case with s > 3/g/2, gas drag is given 
by Stokes' law and the stopping time is given as 

t =l-PL^^ = 1.5x10-^ (-^4^] 7_^^V^_y[year], (9) 
3pg(z)/gC, \lA[cm]J \3[gcm-^]J\2.2[cm]J J' 

at r = 1 AU and z = 0. In this paper, we adopt the criterion s < 3/g/2 for the validity of the 
Epstein regime while there are also previous studies that use s < 9/g/4 (e.g., Youdin & Shu 
2002). 

In this paper, we assumed that radial motion of dust grains can be neglected (Nakagawa 
et al. 1981). The equation of motion of a dust grain in the vertical direction is given by 

-IT = -7 2, (10) 

where t is the time and Vz is the vertical velocity of the dust grain. We assume that gas are 
not affected by dust grains and gas density is always given by equation (5). In this paper, we 
approximate vertical velocity of dust grains by the terminal velocity (Nakagawa et al. 1981; 
Nakagawa et al. 1986). With setting dvz/dt = in equation (10), we obtain a terminal velocity 
of a dust grain as 

Vz{z) = -tstop^K^Z. (11) 

2.4- The mixture of gas and dust grains 

In this paper, we use a single-fluid approximation for the azimuthal motion. The rota- 
tional velocity of a mixed fluid of gas and dust is given as 



Vcj, 



1 MfL^' 



VK- (12) 



In equation (12), p^ is the dust density, vk is the circular Keplerian velocity, and 

IHj^dlnP 
4 r"^ alnr 

where P = Cg^Pg is the gas pressure. Using equations (2), (5), (6) and (7), we have 

^ = Sf-y = l-8xlO-^fTT^7TTV«l- (14) 



16 V r / ■ Vl[AU] 
2.5. Richardson number 

To discuss the possibility of the shear-driven turbulence, we calculate Richardson number 
that is the indicator of KHI (Chandrasekhar 1961). The Richardson number is given by 
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Ri = - 



dz dz 



dz 



(15) 



where is the gravitational acceleration given by = Hk^z (Sekiya & Ishitsu 2001). Equation 
(15) is based on an assumption that the disk is composed of an incompressible one component 
fluid for simplicity. We use equation (15) as an indicator of KHI for simplicity. If Ri is smaller 
than the critical value Ric, KHI is expected to be induced in the protoplanetary disk, and it is 
expected that the laminar flow of the mixed fluid in the disk becomes turbulent. Chandrasekhar 
(1961) showed that Ric = 0.25, but Gomez & Ostriker (2005) and Johansen et al. (2006) showed 
that the inclusion of the Coriolis force yields a much higher critical Richardson numbers. In 
this paper, we use a value Ric = 0.25, but in the subsequent section we will discuss the case 
with Ric = 0.8 (Johansen et al. 2006). 

By seeing equations (12) and (15), Ri depends on the distribution of dust density. In this 
paper, we calculate Richardson number for the dust density given by numerical calculations in 
each time. 

2. 6. The growth and sedimentation of dust grains 

The growth and sedimentation of dust grains is described by 
d d 

— n(m, z) + -Q^ [n(m, z)v^ (m, z)] 

POO 

= —n{m,z) / A{m,m' ,z)n{m' ,z)dm' 



1 Z"™ 

H — / [A{m — m' ,m' ,z)n{m — m' ,z)n{m' ,z)dm'], (16) 

2 JO 

where n{'m,z)dm is the number density of the dust grains with mass between m to m + dm at 
the height z and A{'m,m' , z) is the coalescence rate for two dust grains with m and m' at z. 
The symbol n{m,z)dz gives a mass function in z to z + dz. The dust density at z, Pd{z), is 
given by 

roo 

Pd{z) = / mn{m,z)dm. (17) 
Jo 

As for velocity which induces dust-dust collisions, velocities generated during sedimen- 
tation and the thermal Brownian motion are considered. The relative velocities of two dust 
grains due to sedimentation Avg and to the thermal motion Avb is given by 

Avs = \vz{m,z) -Vz{m',z)\, (18) 



Avj, = /k^J- + ^, (19) 
respectively. We assume that the coalescence rate A{m,m',z) is given by 

A{m,m',z) = 7r{s + s'f{Avs + AvB)ps, (20) 
where Ps is the coalescence probability and Ps is assumed to be 1 as in Nakagawa et al. (1981). 
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For simplicity, we ignore the bouncing barrier (Zsom et al. 2010) and the electrostatic barrier 
(Okuzumi 2009). 

Integrating equation (16) with respect to m, we have 
d f°° d f°° 

— / mn{m,z)dm + — / mn{m,z)vz{'m,z)dm = 0. (21) 
at Jo oz Jo 

Using equations (17) and (21), we have 

g^Pdiz) + —[pd{z)v,iz)]=0, (22) 

which is the continuity equation for dust grains treated as fluid, where the mean sedimentation 
velocity of dust fluid at z is given by 

_ J^Vzjz) mn{m,z)dm 

"^zyz) — — - — - . \^Z6) 

Jq mn[m,z)am 

3. Sedimentation of dust grains without growth 

In this section, in order to clarify the effect of dust growth, we first consider the case at 
r = 1 AU without the growth. 

3.1. The initial condition 

In this paper, dust density pd{z) is assumed to be symmetric to the midplane, so we 
examine only the region of 2; > 0. As initial condition, the initial density of dust is assumed to 
be 



Pd[z = ^ exp 



(24) 



z 

where is time dependent scale height of dust density profile with = Hg at t = 0. 

In this paper, we assume the same initial functional form for Pg{z) and Pd(^); and we 
neglect the dependence of vk on z. Note that Ri = 00 at t = from equation (15), i.e., initial 
state is stable against KHI. This is because the rotational velocity of the mixed fluid of gas and 
dust grains at t = is independent of z (cf. equation (12)), i.e., dv^/dz = in equation (15). 

3.2. Sedimentation of dust grains with single size 

First, we consider the sedimentation of dust grains with a single size. We assume that all 
dust grains are small enough and that the drag force is given by Epstein's law. The characteristic 
time scale for sedimentation of dust grains is given by 

_ z 1 , , 

\Vz{z)\ tstop^K 

From equations (8) and (25), it is shown that the characteristic time scale of sedimentation is 
much larger than the Keplerian period tscd^K ^ 1 in the case when the stopping time is much 
smaller than the Keplerian period tstop^K <^ 1. In such a case with a single size without dust 
growth, profile of dust density evolves in a self-similar manner (Garaud & Lin 2004). Here, 
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Fig. 1. Dust density and Richardson number at r = 1 AU for the case when dust grains have a single size 
with /d = 1- (a) The distribution of dust density (sohd line). The abscissa shows the dust density Pd{z) 
in unit of pdo(O)- The ordinate shows z coordinates in unit of Hg. The critical density pc is also drawn 
(dotted line), (b) The distribution of Richardson number (solid line). The abscissa shows Richardson 
number. The ordinate shows z coordinates in unit of Hg. The critical value Ric is also drawn (dotted 
line). 

during sedimentation of dust, the distribution of gas density assumed to be constant remain as 
equation (5). When the time evolution of dust density proceeds in a self-similar manner, dust 
density profile in t > is also given by equation (24) with temporally decreasing ifd(^) with 
< i/d(^) < Hg. Using equations (5), (12) and (24), we have 



Substituting (14) and (26) into (15), an analytical formula for Richardson number can be 
derived. 

Figure 1 shows results of distribution of dust density and Richardson number when 
Pd(0) ~ 160pdo(0) for the case with /a = 1, where is 0.0062ifg. We define the initial value of 
dust density at the midplane as Pdo(O) = Sd/ {^H^. In Figure 1, it is seen that dust density at 
the midplane is much smaller than the critical density pc = 6.3 x lO^pdo(O), and that Richardson 
number is smaller than Ric = 0.25 around the midplane {z/H^ ^ 0.0075). Thus, KHI is expected 
before GI in this case with /g = /d = 1 at r = 1 AU. 

Assuming equilibrium condition for KHI, previous studies have shown that GI tends to 
occur if dust abundance in the protoplanetary disk is much larger than MMSN model (Sekiya 
1998; Sekiya & Ishitsu 2001). Now using non-equilibrium time dependent density profile during 
sedimentation, we consider the possibility for GI in the case with large /d. We seek the condition 
of dust abundance /d by which GI occurs before the onset of KHI. In Figure 1 (a), it is seen that 
the distribution of dust density has a maximum value at the midplane. When /d is larger than 
1, the dotted line in Figure 1 (a) moves to left because the abscissa is inversely proportional 
to Pdo(O) oc /d. Dust density at the midplane is the first to reach the critical density for GI 




(26) 
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because time development of dust density proceeds in a self-similar manner. A characteristic 
Ha when dust density at midplane attains pc can be obtained from equation (24). We define 
the characteristic as He- With substituting 0, and pc for z, Hg and pd(-z) in equation 
(24), respectively, He is given by 

He = = 1.6 X IQ-'UHg (X /d. (27) 

From Figure 1 (b), it is found that the distribution of Richardson number has a local minimum 
value at z/Hg = 3.5 x 10~^. The height where distribution of Richardson number takes the local 
minimum value is defined as Zc- At z = z^, from dR.i/dz\z=z^ = 0, we find 

Pg{z,) = 2pA{ze), (28) 

with assuming H^/H^ ^ 1. From equation (28), Zc is given as 

2Sd Hg 



Zr. 



In 



Hg. (29) 



Eg Hd^ 

Richardson number a.t z = z^ is given by 

with assuming H^/ Hg<t^l. From equations (27), (28) and (30), the condition of dust abundance 
/d that is necessary for GI to occur before KHI is derived as 



/d 



8vr^ , 
—Ri(z 
27 ^ 



r]rpe > 6.6 x 10^ (31) 



at r = 1 AU and /g = 1. Note that Ed oc fd- 

In Figure 2, results for the case with fd = 6.6 x 10^ are shown. In Figure 2, it is seen that 
dust density at the midplane indeed attains the critical density pc and that Richardson number 
remains marginally larger than the critical value Ric- Thus, in the case when dust grains have 
a single size and don't grow at r = 1 AU with /g = 1, GI is expected to occur before KHI only 
if the dust surface density is much larger than the gas surface density. 

This tendency is similar to the result of Sekiya (1998), but the value of fd in this paper 
(/d = 6.6 X 10^) is larger than the value in Sekiya (1998) (/d = 16.8). The value of fd in Sekiya 
(1998) corresponds to the condition that the protoplanetary disk is in a quasi-equilibrium state 
for KHI. On the other hand, the value in this paper corresponds to the condition that GI occurs 
before KHI during sedimentation and that is more stringent. Therefore, the value of fd in this 
paper is larger than the value in Sekiya (1998). 

3. 3. The effect of size distribution 

Next, we consider the sedimentation of dust grains with initial size distribution but 
without growth. The result given by equation (31) is independent of dust radius as long as the 
stopping time is given by Epstein's law. Now we consider the effect of initial size distribution of 
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Fig. 2. Same as Figure 1, but for the case with /d = 6.6 x 10^. 

dust grains. Characteristic time scale for sedimentation of dust grains in equation (25) depends 
on size of dust grains. Thus, even in the same time, scale height of dust density profile for 
different size of dust grains is different. Therefore, it is expected that total dust density profile 
will change from initial Gaussian profile in the case when dust grains have a size distribution. 
We consider kinds of dust grains with different radius. For simplicity, we consider the case 
when s is given by integer multiples of the minimum radius of dust grains Sq and when the 
maximum value Smax is given by iVdSo. 

In order to derive the scale height H^^s of dust density profile for dust grains with a 
radius s, we assume that Pg{z) is uniform for simplicity. The stopping time of dust grains 
^stop = (Ps/Pg)(s/cs) is proportional to s and is constant with z. The vertical velocity of dust 
grains Vz{z) oc tstop^; is proportional to s and to z. From the time evolution of z{t) with different 
radius of dust grains with z = Hg at t = 0, the formula for H^.s can be derived as 



TT / IT \ */*0 

Hg \ Hg j 



(32) 



In equation (32), it is seen that H^ s is determined only by i^d.so ^"^^ ^/^o instead of s because 
radii of all dust grains are normalized by the smallest dust grains. 
The distribution of total dust density pd{z) is given as 

Pci(^) = E Pd(s,^), (33) 

s/so=l 

where pd{s,z) is the density of dust grains with radius s at z. For simplicity, initial size 
distribution is assumed to be a power law of dust radius, 

nd{s,z) =nd{so,z)(^y^ , (34) 

where 71^(3, z) is the initial condition of the number density for radius s at z with a power index 
p. We assume p = —3 for simplicity. In the case with p = —3, s^nd{s,z) is equal to So^nd{so,z) 
from equation (34). Then, initial condition of pd{s,z) is given by 
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Fig. 3. Dust density and Richardson number at r = 1 AU for the case when dust grains have a size 
distribution with /d = 50. Lines, abscissas and ordinates show the same as ones of Figure 1. 

4 4 

Pd{s,z) = -7ipsS^nd{s,z) = -7rpsSo^nd{so,z) = pd{so,z). (35) 

In the case with p = —3, equation (35) shows that initial density of dust with different radius 
is the same. In this case, pd{s,z) is given by 

In the case without dust growth, pd{s^z) with different s evolves independently, with different 
ifd,s- From equations (32), (33) and (36), we can derive analytical solutions of pd{,z). We 
assume Nd = 1000. 

Figure 3 shows dust density and Richardson number for the case with /d = 50. In Figure 
3, it is seen that dust density at the midplane attains the critical density pc, and that Richardson 
number remains larger than the critical value. This demonstrates that dust fraction 50 times 
larger than MMSN model induces GI before KHI in the case without dust growth. Note that 
this dust fraction is about 10 times smaller than that in Figure 2. We checked the dependence 
on A'^d and found that this is the case with for Nd ^ 10. Thus, it can be suggested that KHI 
tends to be inhibited before GI if dust grains have size distribution even with the same fd- The 
reason is that the vertical gradient of dust density becomes smaller as a result of the continuous 
size distribution. 

Figure 4 shows dust density in the midplane at the onset of KHI for the both cases 
without and with size distribution. The symbol Pd,KH is the dust density in the midplane at the 
onset of KHI. It is seen that dust density in the midplane at the onset of KHI increases with 
increasing the dust abundance in both cases with and without size distribution. Especially, it 
is seen that the condition pd,KH = Pc is attained by smaller dust abundance fd in the case with 
size distribution fd > 50 than without size distribution fd ^ 700. 

Above results are based on Ri^ = 0.25 (Chandrasekhar 1961). However, Johansen et 
al. (2006) showed that Ric = 0.8. We also investigated the case with Ric = 0.8. In the case 
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Fig. 4. The dust density in the midplane and in r = 1 AU at the onset of KHI for the case without size 
distribution of dust grains (triangles) and for the case with size distribution of dust grains (circles). The 
abscissa shows dust abundance /d, and the ordinate shows the dust density pd,KH in unit of p^. The 
approximated curves are also drawn (dotted lines). 

without size distribution, the required abundance of dust for GI is about 1.2 x 10^. On the 
other hand, in the case with size distribution, the required abundance is smaller than 100. This 
abundance is much smaller than that for the case without size distribution, and this result is 
qualitatively same as that for the case with Ric = 0.25. Thus, below, we concentrate on the 
case with Ric = 0.25. 

In the case without dust growth, results are summarized as follows: 

1. In the case when the abundance of dust grains is given as MMSN model, KHI is expected 
to occur when dust density at the midplane is still much smaller than the critical density 
for GI. 

2. GI tends to occur if the abundance of dust grains is larger. 

3. If dust grains have the initial size distribution, the required abundance of dust for GI has 
the possibility to be smaller than that in the case without size distribution. 

Above results are based on the assumption that dust grains don't grow and that initial 
size distribution is assumed. However, in the actual protoplanetary disks, dust grains are ex- 
pected to collide mutually and to grow. If dust grains grow, the size distribution and the largest 
radius of dust grains will change with time. In order to understand the actual condition of dust 
abundance for GI, size distribution which is consistent with dust growth during sedimentation 
is desirable. Next, we consider the case with dust growth. 

4. The case with dust growth 

4.1. The initial condition and the numerical method 

We investigate sedimentation and growth of dust grains at r = 1 AU in the case with 
/g = 1. We solve equation (16) numerically. We assume that the internal density of dust 
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Fig. 5. Results of the test of our numerical method for growth of dust grains. Numerical results (circles) 
are compared with analytic solutions by Trubnikov (1971) (solid lines). 

grains ps is 3.0 gcm~^. As initial condition, we assume that all dust grains have the initial size 
So = 1.0 X 10~^ cm and initial density given by equation (24). 

In numerical calculations, TVD scheme (Roe 1986) is applied to calculate sedimenta- 
tion of dust grains and the method of WS89 (Wetherill & Stewart 1989) is applied for dust 
growth. The mass coordinates rrii are logarithmically divided into 400 mass bins. Our numer- 
ical method is tested using the analytical solution (Trubnikov 1971). Figure 5 shows results 
with our numerical method for growth of dust grains. In Figure 5, it is seen that our numerical 
method has satisfactory accuracy to calculate growth of dust grains. The z coordinates Zj are 
logarithmically divided into 106 spaced grids. The thickness of the nearest grid to the midplane 
is 1.4 X 10~^Hg. Using this spatial resolution, scale height at GI is sufficiently resolved. 

4-2. Possibilities of KHI in the early stage 

Figure 6 shows the distribution of dust density and Richardson number at t = 25 year. 
In Figure 6, although the distribution of dust density changes little from the initial state in 
this short period, it is seen that Richardson number becomes small enough in small z especially 
z/Hg ^ 10""* <^ 1. This rapid decline of Ri did not come out in the case without growth of 
dust grains in §3. We suppose that growth of dust grains is the origin of this rapid decline in 
Richardson number at z/Hg <^ 1 and at t/tscd ^ 1- 

Equations (12) and (15) show that Richardson number is a function of the gradient of 
dust density dpd{z)/dz, so we now think about dpd{z)/dz. First, we think about the case 
without growth of dust grains to evince the effect of dust growth. We think about the case 
when all dust grains have a single size as initial condition because we assume that all dust 
grains have an initial radius sq as initial condition for the calculation in the case with dust 
growth. In the case without growth and size distribution of dust grains, the typical radius of 
dust grains s{z), which is defined as 
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Fig. 6. (a) The distribution of dust density at r = 1 AU and at t = 25 year. The line, the abscissa and the 
ordinate show the same as ones of Figure 1 (a), (b) The distribution of Richardson number at t = 25 year. 
The abscissa shows Richardson number. The ordinate shows z coordinates in unit of _ffg. The critical 
value Ric is also drawn (dotted line). 

_ _ smn{m,z)dm 
Jo mn[m,z)am 

is equal to sq and is independent on z. Garaud & Lin (2004) shows that time dependence of 
dust density proceeds in the self-similar manner in the case without growth of dust grains. In 
this case, the gradient of dust density is given by 

'-^-^--^ (-) 

In order to compare with the case with growth of dust grains, we think about the sedimentation 
velocity for the gravity center of the dust system at z, Vz{z), which is derived from equations 
(11), (23) and (37). At z/H^ < 1 and at t/tscd < 1, Vz{z) is given by 

v.iz)^-^'-^n^z. (39) 

In the case without growth of dust grains, Vz{z) oc 2; at z/Hg <^ 1 because s{z) = sq. 

Second, we think about the case with growth of dust grains. In this case, the typical 
radius of dust grains s{z) is the linear function of z at z/Hg ^ 1 and at t/tscd ^ 1 owing to 
collisions due to sedimentation (see Appendix for the reason). From equation (39), Vz{z) has 
a second-order term of z because the typical radius of dust grains is the linear function of z. 
From comparing the functional form of Vz{z) with that in the case without dust growth, at 
z/Hg<^l and at t/tsed <^ 1, it is suggested that the gradient of dust density is given by 

^^ = 5,z + 5,^5^ (z-^O), (40) 

where 5i(< 0) and 52(> 0) are appropriate values. 

Figure 7 shows the distribution of the gradient of the dust density at z/Hg ^ 1 and at 
t = 25 year for cases with and without dust growth. From Figure 7, it is confirmed that the 
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Fig. 7. The distribution of the gradient of dust density in r = 1 AU and z/H^ <C 1 at i = 25 year. 
The abscissa shows the gradient of dust density dpd{z)/dz in unit of pdo(0)/ffg- The ordinate shows z 
coordinates in unit of Hg. Circles show the case with dust growth and triangles show the one without 
dust growth. 

distribution of the gradient of dust density is approximated by equation (40) with ^2 > in the 
case with dust growth. From Figure 7, we can approximate dpd{z)/dz ~ [—2.0{z/Hg) + 1.0 x 
10~^][pdo{0) / Hg]. This indicates that the distribution of dust density has a local maximum 
value at z/Hg ~ 5 x 10~^ in the case with growth of dust grains. This can be confirmed by 
seeing Figure 8 that shows the distribution of dust density in z/ Hg <^ 1 and at t = 25 year in the 
case with dust growth. By approximating Pg{z) as Pg(0) and Pd{z) as Pdo(O), the Richardson 
number for z/Hg<^l is given by 

•16\V r \ VSA"^ /S, 



z 



2Sg z dpd{z) Hg 



2z ^ dpdjz) Hg "' ^ 
^^g dz pdo(O) 



(41) 

Ed Hg dz pdo(0)J L^g dz pdo(0)J ^ ^ 

Figure 9 shows the distribution of Richardson number at z/Hg '^1 and at t = 25 year in the 
case with growth of dust grains. In Figure 9, it is seen that numerical solutions of Richardson 
number is sufficiently close to the approximate equation (41). In Figure 9, it is seen that 
Richardson number drawn by the solid line is smaller than the critical value Ric around the 
midplane at t = 25 year. 

However, it is doubtful whether KHI occurs. In Figure 8, it is seen that the distribution 
of dust density has a local maximum value at z ^ 0. Under the local maximum point, the 
Rayleigh- Taylor instability (RTI) is suspected to occur because the distribution of the gradient 
of dust density becomes positive (Chandrasekhar 1961). If RTI occurs, the distribution of dust 
density in this region is expected to be adjusted as to be constant (Sekiya & Ishitsu 2001). If 
it is assumed that the growth rate of RTI is larger than that of KHI, and that the distribution 
of dust density becomes constant around the midplane, the gradient of dust density becomes 
zero near the midplane. In this case with assuming that gas density is given by equation (5), 
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Fig. 8. The distribution of dust density in r = 1 AU and z/ ^ 1 at t = 25 year in the case with growth 
of dust grains. The Una, the abscissa and the ordinate show the same as ones of Figure 1 (a). 
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Fig. 9. The distribution of Richardson number in r = 1 AU and z/Hg <C 1 at t = 25 year in the case with 
growth of dust grains. Circles, the dotted line, the abscissa and the ordinate show the same as ones of 
Figure 6 (b). The solid line is given by equation (41). 

Richardson number near the midplane is given by 

This indicates that KHI does not occur. Thus, below, the possibihty of KHI near the midplane 
in the early phase as indicated in Figure 6 (b) is not considered further. 

4-3. The dust density at the onset of KHI for the case with dust growth 

Figure 10 shows dust density distribution at the onset of KHI in the case with dust 
growth with /d = 1. In Figure 10 (a), distribution of dust density has a local maximum value at 
z = zkt ~ 10~'^-ffg, and there is the region where the gradient of dust density becomes positive. 
As discussed in §4.2, in this case, the distribution of dust density is expected to be adjusted to 
be constant by RTI in the region z ^ zbt- Assuming this RTI, we modify the distribution of 
dust density as Figure 10 (b) with mass conservation. Note that this treatment for RTI is crude 
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Pd(z) / Pc Pd(z) ' Po 

Fig. 10. Dust density in r = 1 AU at the onset of KHI for the case with dust growth and /d = 1. (a) The 
dust density obtained by numerical calculations, (b) The dust density adjusted as to be constant (solid 
line) and the dust density not adjusted (dotted line). The dust density not adjusted is equal to the one 
obtained by numerical calculations. For both of Figure (a) and (b), the abscissas show the dust density 
Pdiz) in unit of a-nd the ordinates show z coordinates in unit of iJg. 

and more accurate treatment should be addressed in the future. Hereafter, same modification 
is always applied for the density near the midplane. 

Figure 11 shows the modified dust density profile at the onset of KHI. Figure 11 cor- 
responds to Figure 4 with dust growth. Since the case with large fd is numerically expensive, 
results only for the case with /d < 4 are plotted. In Figure 11, it is seen that dust density at 
KHI, Pd,KH; increases with increasing dust abundance for the case with < 2. On the other 
hand, in the case with /d > 2, it is seen that dust density at KHI decreases with increasing dust 
abundance. This tendency is qualitatively different from Figure 4. In Figure 11, our results 
show that dust density in the midplane at the onset of KHI for the case with /d = 4 is about 
the same as that for the case with /d = 1. As the physical origin of the decline of pd,KH for 
/d > 2, we consider the difference of property of gas drag. After dust grains grow, the law of 
gas drag changes from Epstein's law to Stokes' law. Figure 12 shows the mass function at the 
onset of KHI at the maximum density for the case with /d = 1 and /d = 4, respectively. For the 
case with /d = 1, the mass function has the peak at m/mo ~ 10^^ and the typical size of dust 
grains is 0.9 cm. On the other hand, for the case with /d = 4, the mass function has the peak 
at m/rrio ~ 10^^ and the typical radius of dust grains is 4 cm. In the case with s > 3/g/2 = 2.2 
cm, the stopping time is given by Stokes' law. By comparison of these results, it is suggested 
that the decrease of dust density at KHI for /d > 2 in Figure 11 originates from the change of 
the law of gas drag due to dust growth. 

To confirm this possibility, for reference, we recalculate the evolution using Epstein's 
law for all size. Solid line in Figure 13 shows the dust density at the onset of KHI for this case. 
It is clearly seen that the dust density at the onset of KHI increases with increasing /d. By 
comparing two lines in Figure 13, it is clear that the physical origin for the decline of Pd,KH for 
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Fig. 11. The modified dust density in r = 1 AU at the onset of KHI for the case with dust growth. The 
abscissa and the ordinate show the same as ones of Figure 4. 
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Fig. 12. (a) The mass function in the case with /d = 1 at z = 1.5 x lO^'^iJg and r = 1 AU where the 
distributions of dust density takes a local maximum value at the onset of KHI. (b) The mass function in 
the case with /d = 4 at z = 3.8 x 10~^Hg and r = 1 AU where the distributions of dust density takes a 
local maximum value at the onset of KHI. In both of diagrams (a) and (b), The abscissa shows the dust 
mass m in unit of niQ. The ordinate shows the mass function mn(m,z) in unit of Pdo(0)/TOo- Dotted lines 
show the border line, where the stopping time changes from Epstein's law to Stokes' law. 

/d > 2 is the change of gas drag from Epstein's law to Stokes' law. Therefore, it is significant 
for us to take into account the dust size dependence of the stopping time as well as dust growth 
when we investigate the shear-driven turbulence in the protoplanetary disk. 

5. Discussion 

5.1. Dependence on the heliocentric distance 

We estimate the dependence of the required abundance of dust for GI on the heliocentric 
distance r. First, we consider the case without dust growth. For simplicity, we consider the 
case when all dust grains have a single size. The scale height of dust density profile at the 
onset of KHI can be obtained from equation (30), and we define the scale height as i^KHi- With 
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Fig. 13. The modified dust density in r = 1 AU at the onset of KHI for the case with dust growth in 
cases of using both Epstein's law and Stokes' law (dotted line and squares), and using only Epstein's law 
for all size (solid line and crosses). The abscissa and the ordinate show the same as ones of Figure 4. 

substituting Ric and i^KHi for R-K^ = ^c) and in equation (30), respectively, i^KHi is given 
by 

Hk„,^(^)*„-^1.0x10-//.(^)\ (43) 

Equation (43) shows that -^khi is independent on /d. The dependence of H^. (defined in equation 
(27)) on r is given by = 1.6 x 10-Vd6cc^rg(r/1 [AU])^/^ and we have 

:^ = 6.6xl02(^^^"'(^^^"\ (44) 

The parameter of condensed water ice ^icc is given by 1 (at r < 2.7 AU) or 4.2 (at r > 2.7 
AU) (Hayashi 1981; Hayashi et al. 1985). Equation (44) shows the possibihty that the required 
abundance of dust for GI is smaller outside the snow line than inside. Next, we consider the case 
with dust growth. At the onset of KHI, most of the dust grains with the typical size calculated 
in §4.3 (defined in equation (37)) have settled near the midplane. Thus, the typical size of dust 
at the onset of KHI can be obtained as the size of dust grains which have settled from high 
altitudes to the midplane. As Safronov (1969), we have dm = Airs^psds ~ —ps'^{2sYpA{z)dz. 
Assuming that s = sq sX z = +oo and that s = Sf ^ sq at 2; = 0, we derive Sf as 

The typical sizes of dust at the onset of KHI calculated in §4.3 were 0.9 cm and 4 cm in the 
case with /d = 1 and /d = 4, respectively, and consistent with equation (45). Equation (45) 
shows that the typical size of dust at the onset of KHI is smaller for larger r. On the other 
hand, the mean free path of gas molecules scales as /g oc pg"^ oc ifgSg"^ oc r^^/^, and /g is larger 
for larger r. Thus, the gas drag force tends to be given by Epstein's law in the outer region of 
protoplanetary disk even for a large /d- The Solid line in Figure 13 shows that the dust density 
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at the onset of KHI is larger for larger /d in the Epstein regime. Therefore, the outer region of 
the protoplanetary disk might be more suitable for GI than the inner region. 

Above discussion is different from that in Takeuchi et al. (2012) that showed that the 
inner part of the protoplanetary disk is more suitable for GI. Takeuchi et al. (2012) is based on 
the different condition from ours. Takeuchi et al. (2012) used the scale height of the dust layer 
which is determined by the balance between sedimentation and diffusion of dust grains, while 
we used the scale height of dust density profile at the onset of KHI. 

5.2. The linear analysis of the shear instability 

In this paper, we used Richardson number as an indicator of KHI in order to discuss the 
possibility of the shear-driven turbulence. However, although KHI occurs, there is a possibility 
that some dust grains continue to settle toward the midplane if the shear-induced turbulence 
is weak. An analysis with only Richardson number is insufficient in order to understand the 
effect of the shear-induced instability in the protoplanetary disk. Sekiya & Ishitsu (2001) 
and Michikoshi & Inutsuka (2006) calculated the growth rate of the shear-induced instability 
by solving the linear perturbation equations. In these studies, the growth rate of the shear- 
induced instability depends on the assumed distribution of dust density which is not consistent 
with the formation process of dust layer. Thus, it would be better to calculate the growth rate 
of KHI as well as RTI with the consistent density distribution that is given in this paper. This 
will be addressed in the forthcoming paper. 

5.3. Possibilities of the streaming instability and the fractal growth of dust 

Youdin & Goodman (2005) and Johansen & Youdin (2007) showed that the dynamics 
in the midplane is dominated by the streaming instability. Bai & Stone (2010) showed that 
dust grains with = f^K^top > 0.01 trigger the streaming instability before KHI. In our model, 
at r = 1 AU, Ts > 0.01 corresponds to s > 2 cm. For our calculations, for the case with dust 
growth and /d > 2, the typical size of dust grains at the maximum density is larger than 2 
cm. Therefore, in the dust layer governed by Stokes' law, the streaming instability would occur 
before KHI and GI. 

In this paper, we assumed compact and spherical dust grains. However, since dust 
grains grow due to dust-dust collisions, large dust grains are aggregates of small dust grains. 
Both laboratory and numerical experiments show that aggregates are not at all compact and 
spherical, but have a fluffy structure (Wurm & Blum 1998; Kempf et al. 1999; Okuzumi et al. 
2009). With the same mass, the radius of the fractal aggregate is larger than the radius of the 
compact and spherical dust grain, and varies with the porosity of the aggregate. Therefore, it 
is essential to investigate the evolution of porosity during the growth of dust aggregate. 
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6. Conclusion 



In this paper, we considered the sedimentation and the growth of dust grains, and we 
discussed the possibility of the shear-driven turbulence during dust sedimentation with and 
without dust growth. We assumed that the gas component is in hydrostatic equilibrium and 
that gas are not affected by motion of dust grains. Dust grains are assumed to be compact and 
spherical, and radial motion of dust grains are neglected. We used a single-fluid approximation 
for the azimuthal motion. For dust-dust collisions, it is assumed that dust grains collide and 
coalesce by sedimentation or the thermal Brownian motion. 

Our study shows the following results: (1) The shear-driven turbulence is expected to 
occur before the onset of the gravitational instability in MMSN model at r = 1 AU. (2) In the 
case without dust growth, at the onset of KHl, dust density in the midplane is large for larger 
dust abundance. This tendency holds for the same initial size distribution of dust grains. (3) 
In the case with dust growth, dust density at the onset of KHI decreases with increasing dust 
abundance with f^> 2. This result is qualitatively different from the one in the case without 
dust growth. The reason is that gas drag changes from Epstein's law to Stokes' law for larger 
dust grains which grow up in advance. Thus, we stress that, for the study of shear-driven 
turbulence, the change of the law of gas drag from Epstein's law to Stokes' law as well as dust 
growth is required to be taken into account. 

For the formation of planetesimals, it has been suggested that in order to occur GI 
before inward drift of dust, dust grains have to settle toward the midplane, and it is suggested 
to possible with large dust abundance. However, in this paper, it is suggested that the shear- 
driven turbulence is certain to happen even if dust abundance in the protoplanetary disk is 
larger than MMSN model. In future work, we should develop a more realistic model with 
calculating the effect of KHI and RTI directly. 

We thank Fumio Takahara for fruitful discussion and continuous encouragement. We also 
acknowledge discussion with Sugawara in the early stage of this work. 

Appendix. The typical radius of dust grains in the case with dust growth 

We explain the reason that the typical radius of dust grains s{z) is the linear function 
of z at z/Hg ^ 1 and at t/tgcd ^ 1 owing to collisions due to sedimentation in the case with 
growth of dust grains. We now estimate the mean collision time in the same method used in 
Nakagawa et al. (1981). The mean collision time is given by 

icon = (Al) 

where rid is the number density of dust grains, a is the collisional cross section and Av is the 
relative velocity of the dust-dust collision. We assume that radii of dust grains are given by 
the typical radius and that masses of dust grains are given by the typical mass of dust grains. 
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The typical mass is given by 

m{z) = ^np,[s{z)]\ (A2) 

and we regard that a = vrs^. We treat rid as pd/fn or (Sd/Eg)[pg/m] at z/Hg ^ 1 and at 
t/tscd ^ 1- For the colhsion due to sedimentation, we simply put |s — s'| = s and Av = Avs- 
The mean collision time for sedimentation is defined as tcou.s- At z/Hg<^l and at t/tscd ^ 1, 
tcou.s is obtained by 

In equation (A3), it is seen that tcoii.s is independent of s and tcou.s oc z~^ at z/Hg <^ 1 and 
at t/tscd ^ 1- For the collision due to the thermal motion, we simply put m = m' = rh and 
Av = Avb- The mean collision time for the thermal motion, tcoii,B; is obtained by 

\ 2' 



coll, B 




28(-)'[year], (A4) 
So 



at z/Hg<t^l and at t/tsod <^ 1. In equation (A4), it is seen that tcou.e is independent of z and 
tcoii,B oc s^/^ at z/Hg<^l and at t/tscd ^ 1- 

From equations (A3) and (A4), it is found that tcoii,B < tcoii,s while radii of dust grains are 
relatively small. Therefore, it is suggested that collisions due to the thermal motion is dominant 
as long as radii of dust grains are relatively small. After dust grains have grown, it is expected 
that tcoii,B > tcoii,s, and that collisions due to sedimentation are dominant. The growing speed of 
dust grains in sedimentation is expected to be proportional to z because tcoii,s z~^ . Thus, it is 
supposed that the typical mass of dust grains at z, rh{z), is given by rh{z) = [ai{z/Hg) +a2]mo, 
where ai and 02 are appropriate values. Then, the typical radius of dust grains at z is given 
by s{z) = [ai{z/Hg) + a2]^^^SQ. If aiz/a2Hg <^ 1, it is supposed that s{z) is approximated by 
[a^lz/Hg) + 04)^0 with Taylor expansion. Symbols 03 and 04 are appropriate values. 

Figure 14 shows the distribution of the typical mass and radius of dust grains ai z/Hg^l 
and at t = 25 year. From Figure 14, it is confirmed that s{z) = [a^{z/H^ +a^SQ at z/Hg <^ 1 
and at t/tsed ^ 1- 

Above discussions assume that collisions due to sedimentation become more dominant 
than those due to the thermal motion immediately. However, it is not confirmed that collisions 
due to sedimentation is dominant to growth of dust grains before t = 25 year. We now should 
confirm that this assumption is appropriate for the case that we investigate. 

At z/Hg <^ 1 and t ~ 0, it is expected that collisions due to the thermal motion is 
dominant and tcoii,B is independent of z. In this case, it is supposed that the typical radius of 
dust grains is independent of z, i.e., as in the formula for s{z) is temporally constant. However, 
in a certain time, it is expected that the dominant effect in the growth of dust grains changes 
from collisions due to the thermal motion to those due to sedimentation because of dust growth. 
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Fig. 14. (a) The distribution of the typical mass of dust grains in r = 1 AU at i = 25 year. The abscissa, 
where mtypicai means m, shows the typical mass m in unit of rriQ. The ordinate shows z coordinates in 
unit of Hg. (b) The distribution of the typical radius of dust grains at t = 25 year. The abscissa, where 
stypicai mcans s, shows the typical radius s in unit of so. The ordinate shows z coordinates in unit of Hg. 
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Fig. 15. The time development of as (Circles) in r = 1 AU and z/Hg ^ 1 at t/tsed <C 1. The abscissa shows 
time and the ordinate shows 03 that is derived from fitting s{z) into [a^(z / Hg) + a4]sQ. The approximated 
curve is also drawn (dotted line). 

Therefore, the time for this change can be determined by investigating the time development 
of 03 for s{z). 

Figure 15 shows the time development of a^. From Figure 15, is approximated by 
a3 = 6.3x 10-^(t/l [year])2 + 7.0x 10-3(t/l [year]) - 3.0 x lO'^. This shows that 03 > at t > 4 
year, so it is considered that collisions due to sedimentation dominate in growth of dust grains at 
t ^ 4 year. Therefore, we show that the assumption that collisions due to sedimentation become 
more dominant to growth of dust grains than those due to the thermal motion immediately is 
proper to the case that we investigate. 
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